// THIS FILE FINDS THE LAMBDAS (THE MC OF THE MARGINAL PLANT)
// CALCULATES THE AVERAGE MB (TABLE 1) - (INCL. CO-POLLUTANTS)

set more off
cd "${PATH}/data"

use "entsoe.dta", clear

drop eventtime 

// KEEP OBSERVATIONS FOR 2015 and 2016 ONLY	
drop if yofd(daily) > 2016


// Ajustment for fossiloil: TenneT has one plant running at 1 MWH capacity always from March 2016 onwards 
* Do only consider as marginal source if prodcution >1
replace fossiloil = 0 if fossiloil == 1


// SORT GENERATION SOURCES WITH INCREASING COST
order windoffshore windonshore solar otherrenewable marine hydrowaterreservoir hydrorunofriverandpoundage hydropumpedstorage geothermal waste biomass nuclear fossilbrowncoallignite fossilcoalderivedgas fossilhardcoal fossilgas fossilpeat fossiloilshale fossiloil 


// FIND MARGINAL SOURCE 
local ele_type "windoffshore windonshore solar otherrenewable marine hydrowaterreservoir hydrorunofriverandpoundage hydropumpedstorage geothermal waste biomass nuclear fossilbrowncoallignite fossilcoalderivedgas fossilhardcoal fossilgas fossilpeat fossiloilshale fossiloil" 
local types: word count `ele_type'

forvalues n = 1/`types' {
	local source: word `n' of `ele_type'  
	gen marg_`n' = "`source'" if !missing(`source') & `source'!=0
} 

gen marginal = ""
foreach n of numlist 19(-1)1 {
	replace marginal = marg_`n' if !missing(marg_`n') & missing(marginal)
}


// 50 HZ, special case: Adjustment in fossiloil  //
// REPLACE MARGINAL TECHNOLOGY WITH THE SECOND MOST EXPENSIVE IN 50HZ as OIL production linked to refinery (always running and not responding to prices)
replace marginal = "" if TSO == "50Hertz"
foreach n of numlist 18(-1)1 {
	replace marginal = marg_`n' if !missing(marg_`n') & missing(marginal)  & TSO =="50Hertz"
}
*

drop marg_*


// MARGINAL COSTS ($/MWH)
gen marginal_cost = 0 
replace marginal_cost = 5.1 if marginal == "nuclear" // source: ENTSO-e TYNDP 2018 modeling data as in Golub et al. 
replace marginal_cost = 6 if marginal == "biomass" // allow biomass to be marginal: assign positive value higher than nuclear in line with industry reports
replace marginal_cost = 11.3 if marginal == "fossilbrowncoallignite"  // Brown Coal Report: est. cost of extraction in Germany
replace marginal_cost = cost_MWH_coal if marginal == "fossilhardcoal" // Based on international markets
replace marginal_cost = cost_MWH_ng if marginal == "fossilgas" // Based on international markets
replace marginal_cost = cost_MWH_crude_oil if marginal == "fossiloil" // Based on international markets



// MARGINAL EMISSIONS (t/MWH)

gen marginal_emissions = 0 
replace marginal_emissions = emission_MWH_coal if marginal == "fossilhardcoal"
replace marginal_emissions = 1.148 if marginal == "fossilbrowncoallignite" 
replace marginal_emissions = emission_MWH_ng if marginal == "fossilgas" 
replace marginal_emissions = emission_MWH_crude_oil if marginal == "fossiloil"

* NEW: CO-POLLUTANTS (see "EU-emissions-data_localpollution.do") (t/MWh)
ge marginal_emissions_SOX = 0
replace marginal_emissions_SOX = .0004455 if marginal == "fossilhardcoal"
replace marginal_emissions_SOX = .0007348 if marginal == "fossilbrowncoallignite" 

ge marginal_emissions_NOX = 0
replace marginal_emissions_NOX = .0006804  if marginal == "fossilhardcoal"
replace marginal_emissions_NOX = .00085  if marginal == "fossilbrowncoallignite" 
replace marginal_emissions_NOX = .0002644 if marginal == "fossilgas" 
replace marginal_emissions_NOX =  .0003164  if marginal == "fossiloil"


**********

// Step length in line with technologies 
* Generate groups (main steps) in line with categories used in Supply_Demand_All_TSOs_analysis:
egen step_0 = rowtotal(windoffshore windonshore solar otherrenewable marine hydrowaterreservoir hydrorunofriverandpoundage hydropumpedstorage geothermal)
gen step_1 = nuclear
egen step_2 = rowtotal(waste biomass)
gen step_3 = fossilbrowncoallignite
egen step_4 = rowtotal(fossilcoalderivedgas fossilhardcoal)
egen step_5 = rowtotal(fossilgas fossilpeat)
egen step_6 = rowtotal(fossiloilshale fossiloil)
mvencode step_*, mv(0) override // TSOs without these technologies

// Adjustment of oil plant in 50HZ as above
replace step_6 = 0 if TSO == "50Hertz"

**********
	
// COMPUTE EXPECTED VALUE OF DERIVATIVE OF AS WRT R

merge m:1 TSO time_utc using data_15min
keep if _merge == 3
drop _merge

// Need to merge the cluster for load 
ren daily date
ren hourly hour
merge m:1 TSO date hour using "data15min_cluster_TSO.dta", keepusing(kload)
keep if _merge ==3 
drop _merge


// DERIVATIVES OF REGRESSION MODEL WRT SOLAR (1_AS_regressions.do)
	gen dAS_dR = 5.174 +  (-0.00223) * 2 * solar_gen + ///
		( -9.83e-09) * 3 * solar_gen^2  + (0.000380) * load + ///
		 ( -6.76e-08) * load^2 + ( 0.000000189) * 2 * solar_gen * load if kload == 1 & TSO == "50Hertz"
		
	replace dAS_dR =  6.455 +  (-0.000391) * 2 * solar_gen + ///
		(-3.10e-08) * 3 * solar_gen^2  + ( 0.0000577) * load + ///
		 (-6.62e-08) * load^2 + (5.46e-08) * 2 * solar_gen * load if kload == 3 & TSO == "50Hertz"	
		
	replace dAS_dR =   18.44 +   (-0.00560) * 2 * solar_gen + ///
		(0.000000512) * 3 * solar_gen^2  + ( -0.00260) * load + ///
		  0.000000159 * load^2 + (0.000000185) * 2 * solar_gen * load if kload == 2 & TSO == "50Hertz"			  

	replace dAS_dR = 79.60 +  0.00101 * 2 * solar_gen + ///
		(0.000000101) * 3 * solar_gen^2  + (-0.00558) * load + ///
		 0.000000103 * load^2 + (-9.49e-08) * 2 * solar_gen * load if kload == 1 & TSO == "Amprion"
		
	replace dAS_dR =  -20.40 +  (-0.000250) * 2 * solar_gen + ///
		(0.000000442) * 3 * solar_gen^2  + (0.00184) * load + ///
		 -2.90e-08 * load^2 + (-0.000000117) * 2 * solar_gen * load if kload == 2 & TSO == "Amprion"	
		
	replace dAS_dR =  77.84 + 0.00388 * 2 * solar_gen + ///
		(-6.73e-09) * 3 * solar_gen^2  + (-0.00705) * load + ///
		  0.000000156 * load^2 + (-0.000000155) * 2 * solar_gen * load if kload == 3 & TSO == "Amprion"			  

	replace dAS_dR = 380.1 + (-0.000817) * 2 * solar_gen + ///
		(1.43e-08) * 3 * solar_gen^2  + (-0.0347) * load + ///
		  0.000000788 * load^2 + (3.51e-08) * 2 * solar_gen * load if kload == 1 & TSO == "TenneT"
		
	replace dAS_dR = -3.852 + 0.00105 * 2 * solar_gen + ///
		(-5.31e-09) * 3 * solar_gen^2  + (-0.0000909) * load + ///
		 1.98e-08 * load^2 + (-5.87e-08) * 2 * solar_gen * load if kload == 2 & TSO == "TenneT"	
		
	replace dAS_dR = 42.15 + 0.00161 * 2 * solar_gen + ///
		(-1.91e-08) * 3 * solar_gen^2  + (-0.00512) * load + ///
		  0.000000148 * load^2 + (-6.88e-08) * 2 * solar_gen * load if kload == 3 & TSO == "TenneT"			  

	replace dAS_dR = 130.8 +  0.00561 * 2 * solar_gen + ///
		(-1.23e-08) * 3 * solar_gen^2  + (-0.0329) * load + ///
		 0.00000202 * load^2 + (-0.000000616) * 2 * solar_gen * load if kload == 1 & TSO == "TransnetBW"
		
	replace dAS_dR = 40.25 + 0.0160 * 2 * solar_gen + ///
		(0.00000132) * 3 * solar_gen^2  + (-0.0217) * load + ///
		 0.00000274 * load^2 + (-0.00000375) * 2 * solar_gen * load if kload == 2 & TSO == "TransnetBW"	
		
	replace dAS_dR = 105.7 + 0.0130 * 2 * solar_gen + ///
		(-0.000000160) * 3 * solar_gen^2  + (-0.0314) * load + ///
		  0.00000223 * load^2 + (-0.00000150) * 2 * solar_gen * load if kload == 3 & TSO == "TransnetBW"			  

	  
	  
hist dAS_dR, by(TSO) xlabel(, grid)
// graph export "../figures/hist_dAS_dR_by_TSO.pdf", replace	


**********************************
// TABLE 1: MARGINAL BENEFITS

// LIMIT to solar_gen >0 for table with marginal benefits (standard deviations)
// keep if solar >0 
** keep all observations for reallocation exercise (robust)
**********************************


// COMPUTE EXPECTED VALUE OF MARGINAL COST

gen OC = marginal_cost * omega_15min

egen OC_weighted = total(OC), by(TSO)
sum OC_weighted if TSO == "50Hertz" 
sum OC_weighted if TSO == "Amprion" 
sum OC_weighted if TSO == "TenneT" 
sum OC_weighted if TSO == "TransnetBW" 
lab var marginal_cost "marginal cost ({c 0128}/MWh)"
hist marginal_cost, by(TSO) xlabel(, grid) color(blue%20) ///
	graphregion(lwidth(none) ilwidth(none)) plotregion(lwidth(none) ilwidth(none)) scheme(s1mono)
//graph export "../figures/hist_OC_by_TSO.pdf", replace	

//SD for marginal cost
su marginal_cost if TSO == "50Hertz" 
su marginal_cost  if TSO == "Amprion" 
su marginal_cost  if TSO == "TenneT" 
su marginal_cost  if TSO == "TransnetBW" 


// COMPUTE EXPECTED VALUE OF MARGINAL EMISSIONS
gen avoid_emissions = marginal_emissions * omega_15min
egen avoid_emissions_weighted = total(avoid_emissions), by(TSO)
sum avoid_emissions_weighted if TSO == "50Hertz" 
sum avoid_emissions_weighted if TSO == "Amprion" 
sum avoid_emissions_weighted if TSO == "TenneT" 
sum avoid_emissions_weighted if TSO == "TransnetBW" 
hist marginal_emissions, by(TSO) xlabel(, grid)
// graph export "../figures/hist_avoid_emissions_by_TSO.pdf", replace	


// COMPUTE EXPECTED VALUE OF MARGINAL EMISSIONS (Local Pollutants)
gen avoid_emissions_SOX = marginal_emissions_SOX * omega_15min
egen avoid_emissions_SOX_weighted = total(avoid_emissions_SOX), by(TSO)
sum avoid_emissions_SOX_weighted if TSO == "50Hertz" 
sum avoid_emissions_SOX_weighted if TSO == "Amprion" 
sum avoid_emissions_SOX_weighted if TSO == "TenneT" 
sum avoid_emissions_SOX_weighted if TSO == "TransnetBW" 

gen avoid_emissions_NOX = marginal_emissions_NOX * omega_15min
egen avoid_emissions_NOX_weighted = total(avoid_emissions_NOX), by(TSO)
sum avoid_emissions_NOX_weighted if TSO == "50Hertz" 
sum avoid_emissions_NOX_weighted if TSO == "Amprion" 
sum avoid_emissions_NOX_weighted if TSO == "TenneT" 
sum avoid_emissions_NOX_weighted if TSO == "TransnetBW" 


drop if missing(time_utc)
// replace TSO = "50 Hertz" if TSO == "50Hertz"
saveold  "lambda_15min_locpollution", replace

*******************************************************************************
// DEFINE LOCALS WITH COST OF EMISSIONS (IN EUR/t)
local SCC 31.71   // IN EUR/tCO2
local SOX_cost 34051.2 
local NOX_cost 2021.79
*******************************************************************************

tabstat avoid_emissions_weighted avoid_emissions_SOX_weighted avoid_emissions_NOX_weighted, by(TSO) s(mean N)

gen avoided_emissions_cost = `SCC' * avoid_emissions_weighted 

// Add cost of SOX and NOX emissions:
gen avoided_emissions_SOX_cost =  `SOX_cost' * avoid_emissions_SOX_weighted 
gen avoided_emissions_NOX_cost =  `NOX_cost' * avoid_emissions_NOX_weighted 

tabstat avoided_emissions_cost avoided_emissions_SOX_cost avoided_emissions_NOX_cost, by(TSO) s(mean N)


//SD
gen avoided_emissions_cost1 = `SCC' *marginal_emissions 
su avoided_emissions_cost1 if TSO == "50Hertz" 
su avoided_emissions_cost1  if TSO == "Amprion" 
su avoided_emissions_cost1  if TSO == "TenneT" 
su avoided_emissions_cost1  if TSO == "TransnetBW" 

gen avoided_emissions_SOX_cost1 = `SOX_cost' *marginal_emissions_SOX 
su avoided_emissions_SOX_cost1 if TSO == "50Hertz" 
su avoided_emissions_SOX_cost1  if TSO == "Amprion" 
su avoided_emissions_SOX_cost1  if TSO == "TenneT" 
su avoided_emissions_SOX_cost1  if TSO == "TransnetBW" 

gen avoided_emissions_NOX_cost1 = `NOX_cost' *marginal_emissions_NOX
su avoided_emissions_NOX_cost1 if TSO == "50Hertz" 
su avoided_emissions_NOX_cost1  if TSO == "Amprion" 
su avoided_emissions_NOX_cost1  if TSO == "TenneT" 
su avoided_emissions_NOX_cost1  if TSO == "TransnetBW" 

tabstat dAS_dR, by(TSO) s(mean N)
*tabstat dAS_dR if kload==1, by(TSO) s(mean N)
*tabstat dAS_dR if kload==2, by(TSO) s(mean N)
*tabstat dAS_dR if kload==3, by(TSO) s(mean N)

su dAS_dR if TSO == "50Hertz" 
su dAS_dR if TSO == "Amprion" 
su dAS_dR if TSO == "TenneT" 
su dAS_dR if TSO == "TransnetBW" 


// OVERALL GAINS 
gen MB = OC_weighted + `SCC' * avoid_emissions_weighted + `SOX_cost' * avoid_emissions_SOX_weighted + `NOX_cost' * avoid_emissions_NOX_weighted  - dAS_dR
tabstat MB, by(TSO) s(mean N)

// SD for total
gen MB1 = marginal_cost + avoided_emissions_cost1 + avoided_emissions_SOX_cost1 + avoided_emissions_NOX_cost1 - dAS_dR
su MB1 if TSO == "50Hertz" 
su MB1 if TSO == "Amprion" 
su MB1 if TSO == "TenneT" 
su MB1 if TSO == "TransnetBW" 
